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ABSTRACT 

The sine-Gordon equation is a nonlinear partial differential 
equation. It is known that the sine-Gordon has soliton so- 
lutions in the ID and 2D cases, but such solutions are not 
known to exist in the 3D case. Several numerical solutions 
to the ID, 2D, and 3D sine-Gordon equation are presented 
and comments are given on the nature of the solutions. 

1. INTRODUCTION 

Many problems in science and engineering can be expressed 
as a partial differential equation (PDE). Finding accurate 
numerical solutions is often important for these applications, 
but there are many challenges associated with computing 
these solutions. Nonlinear PDEs can cause numerical meth- 
ods to destabilize[2]. Finding accurate solutions can require 
high resolution, large amounts of memory, small time-steps, 
and a lot of computation time. 

Numerical methods that compute solutions for PDEs must 
approximate derivatives. One way to achieve this is with a 
finite difference; one of the simplest such methods is the 
forward difference given by 



u(x + Sx) - 
Sx 



u(x) 



(1) 



[18| Ch. 1] shows that the finite difference approximation 
will converge to the true value at a rate polynomial in Sx; 
that is, as Sx gets smaller, the approximation gets better 
at a polynomial rate. Getting a very accurate approxima- 
tion of the derivative may require Sx to be extremely small, 
which may exhaust memory requirement or take a lot of 
time, especially in the 3D case. |18| Ch. 4] goes on to 
show that spectral methods converge faster than any poly- 
nomial, so we can get very accurate approximations with a 
relatively larger Sx. This can save a lot of time and memory. 
If periodic boundary conditions are used, we use the Fourier 
transform to compute approximate spatial derivatives. 

This paper considers the sine-Gordon equations in 1, 2, 
and 3 dimensions. The ID problems considered here are all 
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suited to being run on one processor. However, as the di- 
mensionality is increased, several issues make using multiple 
processors necessary. Computing multi-dimensional Fourier 
transforms amounts to computing multiple one-dimensional 
Fourier transforms. These one-dimensional transforms are 
independent and can be computed more quickly in parallel. 
Memory requirements also increase dramatically as the di- 
mensionality increases. Problem sizes of interest, especially 
in 3D, may require more memory than one machine can rea- 
sonably have, so using multiple computers may be necessary. 
Using multiple computers also gives more total cache, which 
generally improves performance. 



2. BACKGROUND 



2.1 Sine- Gordon Equation 

The sine-Gordon equation is a nonlinear partial differen- 
tial equation given by 



Utt 



This equation arises in many different applications, includ- 
ing propagation of magnetic flux on Josephson junctions, 
sound propagation in a crystal lattice, and several others 
discussed in [17]. We consider the 1, 2, and 3 dimensional 
cases given by 



utt — U X x 

Utt Uxx Uyy 

Utt U X x Uyy 



sinw (3a) 
sinw (3b) 



sin it 



(3c) 



The one dimensional sine-Gordon equation is exactly in- 
tegrable. This means a large class of exact solutions can 
be found using the inverse scattering transform, including a 
special class of solutions called solitons [9] . Solitons are non- 
linear effects that lead to travelling waves that are localized 
and permanent. Solitons are interesting physically and as a 
way to check numerical codes. For the sine-Gordon equa- 
tion, some of the explicit exact solutions cannot be easily 
evaluated 13 . For these solutions, it may be may be easier 



to compute the solution numerically and can still give great 
insight. The sine-Gordon equation is known to have soliton 
solutions in the 1 and 2 dimensional casesllOl. 



The sine-Gordon equation has a Hamiltonian E given by 



L\ = I ^|Viij 2 dx 



Ep = J 1 — cos udx 
E = Ek + -Es + Ep 



(4a) 

(4b) 

(4c) 
(4d) 



We give each component of the Hamiltonian a name and 
interpretation[l5]. Ek is the kinetic energy, Es is the strain 
energy, and Ep is the potential energy. E is the total energy 
and is constant. These quantities are useful in evaluating the 
accuracy of numerical schemes and the nature of a solution. 

2.2 Discrete Fourier Transform 

The Discrete Fourier Transform (DFT) can be used to 
approximate derivatives. The DFT and inverse DFT are 
given by 



N 



e ' kXj Vj, k 



N/2 



N 



+ 1, 



u b1 = ^ v k e lkXi , j = l,. 



N 
' ~2~ 

, N 



(5) 
(6) 



fc=-JV/2+l 



Conceptually, we can think of the DFT as decomposing 
v[j] into into a set of orthogonal complex exponentials. Each 
Vj represents the magnitude and phase of its corresponding 
complex exponential. 

For smooth functions, the Fourier coefficients converge 
faster than any polynomial |18| . Using the Fast Fourier Trans- 
form (FFT) algorithm, the DFT can be computed in 0(n log n) 
time [6]. 

Computing the derivative in Fourier space is a simple mul- 
tiplication by the wave number k and i. 



fc=~]V/2+l 
N/2 

i- Y, ikvk^i, j = l,...,N (8) 



2?r 



k=-N/2+l 



2.3 Numerical Method 

Following [7] and [5] , we use a leapfrog method to approx- 
imate the second derivative in time with a central difference 
given by: 



tion by the wave number in Fourier space, giving 
O {{ik x f + (iky) 2 + (ik z )) 



~(.^x H~ fay H - ^z) 



+ 2u n ■ 



(11) 
(12) 



where k x ,k y , and k z are the wave numbers for the x, y, and 
z dimensions, respectively. In the ID (and 2D) case, the k z 
(and k y ) wave numbers are not present. 

We note that because the sine function is nonlinear, sin u 7^ 
sin (u) and we must recompute sin u separately from u for 
every time step. 

Putting everything together, we have 



+ (kx + k y + k' z ) 



5t 2 
-sin u" 



(13) 



Moving terms in Equation 13 around to compute u n+1 
gives 



(2u n - u"- 1 ) (kl + k 2 y + k'i) (2u" + u n 



5t 2 



(14) 



3. IMPLEMENTATION 

3.1 Simulation 

The ID, 2D, and 3D codes are all implemented differ- 
ently to illustrate many programming techniques. The ID 
code is relatively straight-forward FORTRAN implementa- 
tion. The 2D cases of interest are small enough to fit on one 
GPU, so the 2D code is implemented in C using Nvidia's 
CUDA. The FFTs and time-stepping scheme are run en- 
tirely on the GPU. We have found that GPU codes can run 
approximately 50 times faster for this type of code[5]. The 
3D problems of interest use grid sizes that are too large to 
easily run on one computer, so we provide a parallel MPI 
implementation based on [3] . The 2decomp library is used to 
provide an efficient abstraction for performing 3D FFTs [l2] . 

3.2 Visualization 

The ID and 2D plots were generated using Python's Mat- 
plotlib 11 . The 3D plots were generated using Visit [4 . 
Visit is useful for these because it can run in parallel. The 
3D codes can generate a large amount of data, so this feature 
is very important. 



where 8t is the time step size and u n is the approximation 
to the function u at the nth time step. 

We approximate the spatial derivative with a spectral 
method. Differentiation in real space becomes multiplica- 



4. RESULTS 
4.1 ID Breathers 

One type of solution for the sine-Gordon equation is called 
a breather. A breather is a nonlinear mode that is localized 
in space and oscillates with time. In [§], Enz shows that 
these breathers can be interpreted as de Broglie waves; that 
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Figure 1: The interaction between two breathers in one dimension. The middle breather is moving to the 
right and the breather on the right oscillates in place. 1(a) shows before any interaction, 1(b) shows as the 
breathers are interacting, and l(c)| shows after the breathers have interacted. |l(d)| shows how the various 
energies change with time. Note that energy oscillates between potential and kinetic and strain energy 
increases during the interaction between the two solitons, but the total energy remains constant. These were 
generated with A?* = 2 12 , <5i = 0.05, a* = 2, c = 0, 0.9. 



is, a particle that is also a wave. In the ID breather 
is given by [To] [19] 



7 = 



i / cos (7 cos a(t — xc)) 

u(x,t) = 4 tan tanu -V —, tt 

cosh (7 sm/j(i — ct)) 



(15) 

where /1 is a parameter that determines the size and fre- 
quency of the pulse and c is the speed of the pulse. 

|Figure i] shows a numerical solution with a stationary and 
moving breather colliding. After colliding, both breathers 
continue on their original trajectories. 

4.2 2D Breathers 

give initial conditions for stationary and moving breathers 



14 



in 2D. 



i(x, y, t) — —4 tan 1 



A 



breather. 

|Figure 2] shows the results of running the simulation for 
the 2D ring solitons. 

4.3 3D Pseudo-Breather 

A 3D pseudo-breather was constructed based on HI. It is 
called a pseudo-breather because, like a breather, it oscil- 
lates in place. The initial conditions are 



i(x, y, z, to) = 4 3 tan 1 



sin ((f)(t) - kx) sech (X(t)(x - ((t))) sech (X(t)y) 
(16) 



tan 



tan 



v/i- 




sin 


A* 




sin 


v/i- 








sin 







(\/l - At 2 io)sech(^x)) 
(s/l - fi 2 t )sech(fiy)) 



(\/l - ^ 2 to)sech(^2:)) 
(17) 



For k = and f (t) constant this gives a stationary breather, 
otherwise the breather is moving. |Figure 2] shows the simula- 
tion results for initial conditions for a stationary and moving 
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Figure 2: 2(a) and 2(b) show a stationary breather at two different times. 2(e) shows the energy for this 
solution. 12(c) and 2(d) show a moving breather. 2(f) shows the energy for the moving breather. These were 
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2 9 , 8t = 0.01. For the stationary breather, A(0) = 0.1, 0(0) = 0.8, C(0) = 0, k = 0. For the 



moving breather, A(0) = 0.1,0(0) = 0.8, C(0) = 0.3, k 



-0.1. 
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Figure 3: A 3D pseudo-breather, shown with semi-transparent isosurfaces 
place, so the isosurfaces move in and out from the center of the pseudo-breather 
pseudo-breather at different points between an oscillation 
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These were generated with N x = N v = N z = 2 s , St = 0.05, fi 
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3(b)| shows the 



3(a) and 



3(c) shows how the various energies change with 
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-5.4. This was run on Flux with 8 



5. FUTURE WORK 

So far, the 3D codes have only been run on relatively 
small machines with small grid sizes. Further work includes 
moving to a larger machine and larger grid sizes, as well as 
trying more complex interactions between a larger number 
of solitons. 
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